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Abstract 


solwnd-1 . 0 is a three-dimensional compressible MHD code written in For- 
tran for studying the solar wind. Time-dependent boundary conditions are 
available. The computational algorithm is based on Flux Corrected Trans- 
port and the code is based on the existing code of Zalesak and Spicer [1]. 
The flow considered is that of shear flow with incoming flow that perturbs 
this base flow. Several test cases corresponding to pressure balanced mag- 
netic structures with velocity shear flow and various inflows including Alfven 
waves are presented. 

Version 1.0 of solwnd considers a rectangular Cartesian geometry. Future 
versions of solwnd will consider a spherical geometry. Some discussion of 
this issue is presented. 



1 Introduction 


For the purposes of studying 3D, fully compressible MHD flows with the 
capacity to capture shocks a new series of codes based on Flux Corrected 
Transport (FCT) termed solwnd have been developed. The present ver- 
sion of the code considers a rectangular geometry while future versions will 
be cast in spherical coordinates and consider a general frustum geometry. 
Time-dependent boundary conditions are available. These code are tar- 
geted towards addressing three effects previously unapproachable in solar 
wind studies: 


1. Time-dependent inflow boundary conditions which can be imposed as 
well as various combinations of outflow and periodic boundary condi- 
tions. 

2. The capacity to capture shocks, i.e. oscillation-free discontinuities 
spread over not more than a few grid points allowing the consideration 
of severely nonlinear cases and high Mach numbers. 

3. The consideration of spherical geometry allowing the study of the fan- 
ning out of the flow with radius. 


solwnd_l . 0 addresses the first two. The next version, solwnd_2 . 0 , will ad- 
dress number 3. The code solves the inviscid compressible MHD equations in 
conservation form. The FCT procedure is well known [2, 3] and has been em- 
ployed previously for MHD [4, 1]. A purely hydrodynamic version has been 
parallelized to run on several machines [5]. The FCT procedure proceeds by 
first forming a low-order and a high-order flux. Then an anti-diffusive flux 
is formed as their difference. This is subsequently (Bux-) limited to prevent 
the formation of spurious extrema. Finally, the solution is updated based 
on the low-order solution and the anti-diffusive flux. 

The magnetic field is calculated from the electric field via Faraday’s Law. 
The cell-averaged magnetic field is cell-centered, while the electric field is 
edge-based. In this way the divergence condition is met to within numerical 
accuracy. The effect of computing the magnetic field via Faraday’s Law is 
that B is computed in planes and variables are packed and unpacked in two 
dimensional arrays. This reduces storage. 
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The code advances the primary cell-averaged variables of density, the three 
components of momenta, the total energy (internal + kinetic -I- magnetic), 
and the three components of magnetic field. The equation of state is taken 
to be ideal gas with a fixed adiabatic index. All other variables (pressure, 
electric field, low-order solution etc.) are considered temporary. 


2 Equations 


The basic equations of MHD in standard notation are[6]: 


dp_ 

dt 


-I- V • pa - 0, 


'du 

at +u ' Vu 


= -Vp + -(JxB), 

C 


dB 

dt 


-cV x E, 


J = t~V x B, 


J = cr(E + u x B). 


The induction equation, 

dB „ , c 2 2 t^ 

— = V x (u x B) + - — V 2 B. 
dt Ana 

For infinite conductivity, a f (5) implies, 

E + u x B = 0. 

The induction equation becomes, 


(1) 

( 2 ) 

(3) 

(4) 

(5) 


( 6 ) 


(7) 
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— ^Vx(uxB). 


( 8 ) 


The last term can be written, 


V x (u x B) = u(V • B) - B(V • u) + (B • V)u - (u • V)B, (9) 


or 


= V • (Bu - uB). 


( 10 ) 


The dyadic terms are interpreted as 


Vab = (a • V)b + (V • a)b. 


The Lorentz force term in the momentum equation can be written, 


-(J x B) = — -p-B x (V x B). (11) 

c 

The identity, 

Jv(B- B) = (B- V)B + B x (V x B), (12) 

gives, 

i(J x B) = -V(f£) + i(B ■ V)B. (13) 

In addition to the above equations there is the energy equation, one form of 
which is [7], 


De_pDp 
p Dt 


(14) 


Taking u- (2) and B- (8) and adding to the energy equation (14) gives the 
equation for total energy, 
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d 1 _£?2 

Qjipe + 2 f™ 2 + 8^) + V ' K P e + 2 pu2 + p ) u + X = °- ( 15 ) 

Thus the equations as solved in the code become, 

+ V • pu = 0, (16) 

+ V - t(p + + PU" - -^ BB ], (17) 

^(pe + ^pu 2 + |^) + V • [(pe + ipu 2 + p)u + -^u x B x B)] = 0, (18) 

at* 

— : + V-[uB-Bu] = 0, (19) 

with the relations, 


V • B = 0 , 


P = 


(21) 


3 Spherical Coordinates 

For solving these equations in a spherical geometry which future versions of 
solwnd will consider, the code must be written with general metric coeffi- 
cients so that any system of orthogonal coordinates are usable. 

If A = {.4 iXi,. 42X2, A3X3}, the divergence operator in curvilinear coordi- 
nates is: 


V- A = 


hi/12/13 


c\ a 

-^-(h^hzAi) + ^-(^1/13^2) + ^-(hi/i 2 43) 


, ( 22 ) 
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where for spherical coordinates the metric coefficients are 


h\ == 1, /12 = r, h$ = r sin 9. 

So that the divergence operator is 


(23) 


V ' A = i|;(r 2 4) -I- -1-^ A ( Si n 04,) + -A-.±A, 


r 2 dr v ' r sin 989 
The divergence of the dyadic terms is: 


rsin# 8<j) 


<t>- 


(24) 


(V ' AB), = ±l(r*A r B r ) + -LjA (s to04B r ) + ~^ B r) 


r sin# 89 

A g B e + AfpB,p 

r ’ 


r sin# 8<f) 


(25) 


(V • AB)* = \^-{r 2 A r B e ) + -4-^-{sm9A e B e ) + -4“ 
r 2 dr r sm 9 89 r sm 6 d<f> 

AgBr cot 9A^B$ 


(A'pBg) 

(26) 


(V • ABU = + 4ji(VW 


; Atj)Bf cot 9 A ( pB q 

r r 


(27) 


4 Flux Corrected Transport 


solwnd uses Flux Corrected Transport (FCT) as a shock capturing scheme. 
In flux conservative form (for ID) [3] the equation, 


dw df 
dt dx 


= 0 , 
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is 

w i +l = < - 1/2 - -Pi-1/2). 

The tw are interpreted as density, momentum, and the energy; i.e. in ID 
w m {p, pu, pE} and / = {pu, pu 2 4- p, puE + pu}. The dependence of F 
on / is determined by the order of the scheme (see next section). The 3D 
calculation is extended analogously. 

The procedure is as follows, 

1. Form a low-order flux : F^_ j_ x , 2 . 

2. Form a high-order flux : F^^. 

3. Form an anti-diffusive flux: 

Ai+1/2 — -Pi+l/2 — -Pi+l/2- 

4. Form a low-order solution: 

w i = w i ~ ^^+1/2 - ^-1/2)- 

5. Limit the anti-diffusive flux to prevent spurious extrema in 

■^i+ 1/2 = Q+l/2-4i+l/2 0 < Cj + 1/2 < 1- 

Determination of the appropriate C is the critical and complex step. 

6. Update the solution: 

w ? +1 = W i ~ _ A?- 1/2)* 

5 Formation of Fluxes 

In spherical coordinates most of the terms can be written in conservation 
form. However some terms arise that are nonconservative and must be 
treated as source terms for the otherwise conservative equations. The latter 
terms can be treated in a straightforward manner by the inclusion of the 
appropriate metric coefficients in the formation of the various fluxes required 
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in the scheme. The source terms are not flux-limited and are added just prior 
to the solution update. Consider for example the momentum equation, the 
flux conservative form of which is, 

^ + V .[(p + g) {ij+ p««-iBB]. 

In spherical coordinates the x-component becomes, 


+ ii r 2 (p +(m 2 + —) + 

dt r 2 Or ^ r 8ir r sin 9 09 


(sin 9pu$u r ) + 


1 0 
rsin# 0(f> 


(pU^Ur) 


1 0 
Anr 2 dr 


(r 2 B 2 T ) 


Anr sin 9 09 


(sin 9BeB r ) — 


0 


Awr sin 9 d(j> 




+ _ B$ + Bl 

V Airr 

The last two terms, proportional to 1 /r, are non-conservative and have no 
analog in Cartesian coordinates and are treated as source terms in spherical 
coordinates. All the other terms in the above equation are conservative 
and are treated as in the Cartesian case in the FCT procedure. In ID for 
instance, 

w? + 1 =w?-^(F i+ 1 / 2 -F i - 1 /2 ), 

represents the update of a quantity w. The consistent formation of fluxes is 
the key first step in the procedure. For a fourth-order scheme [3], 

7 1 

Fijk = Y 2 ( A?* + fi+l,jk) ~~ ^(/i-lJA; fi+2,jk)‘ (28) 

6 Calculation of B 

The electric field is calculated from the the magnetic field via (7). The 
magnetic field in turn is updated from the electric field via Faraday’s Law. 
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The cell-averaged magnetic field (e.g. k)) is located at the face of the 

cell centered at (i, j, k) (see Figure 1), while the electric field is edge-based. 
Thus the divergence condition is met to within numerical accuracy. Also, 
storage is reduced since B is computed in planes and variables are packed 
and unpacked in two-dimensional temporary arrays. 



Figure 1: Calculation of B from E. 


7 Time stepping 


The code presently uses 2nd order Runge-Kutta (modified Euler) scheme 
for time stepping. With this scheme, 


w n+1/2 = w n ~ =^F{w n ), 
W n+l =w n - t\tF{w n+l ! 2 ). 


The scheme has a weak instability along the imaginary axis, but for the 
length of integrations typically considered here it does not manifest itself. 
The advantage of the scheme of course is that only two time levels of storage 
are required. 
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8 Geometry 


The present code considers a rectangular Cartesian grid. However future 
versions of the code will consider a general frustum (Fig. 2) geometry. This 
has the advantage that the expansion of the flow with radial distance from 
the Sun is naturally incorporated. 


R 2 $sin® 



Figure 2: Frustum geometry of solwnd_2 . 0 


9 Initial conditions 

To simulate one set of model problems relevant to the solar wind, we consider 1 
that the base flow is in the ^-direction, with the direction of shear in the z 
direction and periodic conditions in y and z directions. Then the base flow 

1 This scheme was suggested by A. Roberts and M. Goldstein. 
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Bo x (x,y,z) = B 0 + 8B 0 {y,z), 


(29) 


is given by 2 


or 


B 0x (x , y, z) = Bo + <©o f{y, z), 
where f{y,z) can be chosen, as e.g., 

/(y, z) = cos(2'nk y y)cos(2wk z z). 
We introduce the ratio T)b such that, 

Bq = RMSitooffaz)) ' 

IB 

The transverse magnetic field amplitude is given by, 

y Z “ yj_Bo. 

The initial velocity profile is given by, 


(30) 


(31) 


(32) 


(33) 


v _ f Vo 4- SV o tanh(z - L z /A)/z s if z < L z /2 

x\ iV-, i y Vo — 8Vq tanh(z — ZL z /A)/z s otherwise 

z s is the shear layer thickness parameter. 

We also introduce yy such that 


(34) 


5V 0 = yy V 0 , (35) 

2 To indicate what quantities are true constants we show them in roman type. Thus 
Bo, <5BoVo, <5Vo, are numerical constants. Other constants, not so indicated, are parameters 
such as M and Ma- Spatial dependence is indicated on fields as (x,y,z). Thus Bo is a 
constant while Bo(x,y,z) is an initial field. f(x,y,z) or f(y,z) etc. indicate arbitrary 
functions of the variables, and can indicate different functional dependence from relation 
to relation. 
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<5Vo is determined once Vo is specified. 

Initially the transverse velocities are set to zero, 

V y (x,y,z,0) =0, (36) 


and 


V z (x,y,z,0) = 0. 


(37) 


Now setting, 


Po(x,y,z) = — , 


gives the Alfven speed, 


Va = 


Bp 

(471-po) 1 / 2 


(= Bo, by construction). 


(38) 


(39) 


Since the Alfven Mach number is a parameter Vq can now be specified, 


Vo = M a V a . (40) 

The sound speed is determined from the sonic Mach number, M, 

c s = MV o, (41) 


To obtain pressure balance, if 

p(x, y, z, 0) = + — {po + Sp{x , y,z)) + p tot - + ^-po + ptot, 

o7T 7 07T 7 

canceling common terms, 


p{x, y, z, 0) = po + g ~2 (Pi - ®o (y 5 z ) 2 ) ■ (42) 
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Finally note that since both the sonic and Alfven Mach numbers are param- 
eters, the plasma (3 is a derived quantity, 


B= Vl = M^_ 

<4 Ml' 


(43) 


Table 1 lists the user specified parameters and their Fortran names in the 
solwnd code. 


Parameter 

symbol in text 

Fortran name 

Perturbation B x 

®0 

bx_pert 

Initial B ratio, RMS|^ 

f]B 

eta_b 

Initial V ratio, ^ 

Vv 

eta_v 

Initial transverse B ratio, 

Vx 

eta_perp 

sonic Mach number 

M 

fmach 

Alfven Mach number 

M a 

fmacha 

Bo wavenumber in y direction 

ky 

xky 

Bo wavenumber in 2 direction 

k z 

xkz 


Table 1: Initial condition parameters of solwnd code. Also see section 2 
(page 30). 


10 Time-dependent boundary conditions 

The solar wind can be modeled as a spatially evolving shear layer with time- 
dependent inflow. For the case when the inflow is considered to consist of 
Alfven waves, at each time step these conditions axe imposed at the x = 0 
boundary. Namely, 


B y (0,y, z,t) = fo yz sinwAt, 

(44) 

B z (0, y , z, t) = ® yz cos u>At, 

(45) 

V y (0,y,z,t) = -B y {0,y,z,t), 

(46) 
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V z (0,y,z,t) = -B z {0,y,z,t). 


(47) 


Here u>a is the frequency of the incoming Alfven wave and is an input pa- 
rameter for the code. The density and pressure are fixed as, 


p(0,y,z,t) = p(0,y,z,0), 

(48) 

p(0,y,z,t) = p(0,y,z,0). 

(49) 


Since the code evolves density, momenta and energy, the above imposed 
fields are used to calculate these variables at x = 0. 

11 Test results 

We show three test case results. The first case is a two-dimensional purely 
hydrodynamic calculation. The run was made by making nx=324, ny=5 , 
nz=36; the update on the y — momentum was skipped; and the boundary 
condition in the y — direction was such that all planes in y identical. Figure 
3 shows the vortex roll-up from spatially evolving double shear layer. 

The second test case (Figure 4) considers the triply-periodic geometry of 
solwnd_0 . 2. This temporally evolving shear layer is what is typically studied 
with spectral (Fourier) codes. 

The third test case is with the version solwnd_0.9 with initial conditions 
set forth in the previous sections (version 0.9 contains the same physics and 
algorithm of version 1.0). The parameters of the case are given in Table 2. 
The initial velocity profile is shown in Figure 5. Figures 6 and 7 show the 
fluid and magnetic fields in the x — z plane while Figure 8 show these field in 
the y — z plane. Subsequent figures 9 -14 show cuts of velocity and magnetic 
field as arrow plots so that a sense of direction can be discerned. A streaming 
set of shock waves and magnetic structures are seen to be swept through the 
computational domain. A small subsonic and sub-Alfvenic region is found 
to form near the inflow, no doubt corrupting the flow somewhat. However 
the effect of upstream traveling waves appears to be minimal. Figures 15 
and 16 show the max/min values over the entire computational domain of 
the sonic and Alfvenic Mach numbers, density and pressure. 
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Figure 3: Test case 1. Vortex roll-up. Two-dimensional purely hydrody- 
namic calculation. The resolution is 320 x 32. The vertical scale has been 
exaggerated. 
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t>0 


Figure 4: Test case 2. Triply periodic flow with the indicated parameters 
and configuration. The resolution is 64 3 . The upper panel is at time, t=0, 
while the lower panel is at a later time. 
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Figure 5: Test case 3. Initial double shear layer velocity profile for 

solwnd_0.9 with parameters given in Table 2. The resolution is 256 x 64 2 . 
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Figure 6: Test case 3. Images of density, pressure and y— vorticity, u y , in 
the x — z plane at time t — 5.4 (step 2000) for solwncLO . 9 with parameters 
given in table 2. The resolution is 256 x 64 2 . 


Figure 7: Test case 3. Images of magnetic fields, B x , B y and B z in the x — z 
plane at time t = 5.4 (step 2000) for solwnd-0 . 9 with parameters given in 
Table 2, The resolution is 256 x 64 2 . 
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Figure 8: Test case 3. Images of x— vorticity, density and B x , in the y — z 
plane at time t = 5.4 (step 2000) for solwnd_0 . 9 with parameters given in 
Table 2. The resolution is 256 x 64 2 . 
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Figure 9: Test case 3. Magnetic field vectors in the x — z plane. Left panel 
has components (B x — Bq,B z ). Right panel has components (B x — Bo, 0). 
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Figure 10: Test case 3. Magnetic field vectors in the x — z plane. Left panel 
has components (0, B z ). Right panel has components (0, B y ). 
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Figure 11: Test case 3. Velocity field vectors in the x — z plane. Left panel 
has components (V x — Vo, V z ). Right panel has components (V x — Vo, 0). 
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Figure 12: Test case 3. Velocity field vectors in the x — z plane. Left panel 
has components (0,V Z ). Right panel has components (0, V y ). 
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Figure 13: Test case 3. Velocity field vectors in the y — z plane. Top left 
panel has components (V^, V z ). Top right panel has components (V y , 0). 
Bottom panel has components (0, V z ). 
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Figure 14: Test case 3. Magnetic field vectors in the y — z plane. Top left 
panel has components (B y ,B z ). Top right panel has components (B y , 0). 
Bottom panel has components (0 ,B Z ). 









tl-ne 


Figure 15: Test case 3. Top: Maximum and minimum sonic Mach number 
versus time, Bottom: Maximum and minimum Alfvenic Mach number versus 
time, for solwnd_0 . 9 with parameters given in table 2. The resolution is 
256 x 64 2 . 
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Parameter 

value 


0.1 

» 7b 

0.1 

V± 

0.3 

Vv 

0.3 

M 

1.8 

M A 

1.3 

kx 

0 

ky 

1 

kz 

1 

7 

1.67 

Lx 

4 

Ly 

1 

L z 

1 

Zs 

1 

U A 

13 

A 

10 

B o 

0.5 

Vb 

0.65 

x o 

0.195 

PO 

0.08 

V A 

0.5 

C S 

0.36 

p 

1.9 


Table 2: Initial condition parameters of solwnd code for test case 3. Derived 
quantities are shown below the line. 


12 Version history 


This section records the history of the various versions of solwnd . 


0.0 Original FCT MHD code of Zalesak and Spicer. Non-working. 

0.1 (5-95) Included boundary conditions on cell-averaged B. Allows code 
to work. 

0.2 (5-95) Triply periodic code. Initial conditions introduced: 2 shear lay- 
ers in x— direction with tanh profile. 64 3 grid test cases computed 


28 




with supersonic and super- Alfvenic initial conditions. 

0.3 (6-95) Inflow-outflow in ^-direction and outflow in y, z-directions code 
produced. Because the code computes current from B— field, and then 
updates current to update B field, the n x — 3 point does not get a 
definition for B. Hence B made continuous at n x — 3. Similarly the 
end planes in the x — y and x — z planes need to be made continuous. 

0.4 (6-95) Inflow-outflow in ^-direction and periodic in y, ^-directions code. 

0.5 (7-95) Time-dependent boundary conditions added. Initial conditions 
corresponding to a pressure balanced double shear layer with sinu- 
soidal transverse B and incoming Alfven waves added. 64 3 and 96 3 
computations done. First results of spatially evolving cross-helicity 
shown. Presented in Goddard report ????. 

0.6 (1-96) Modified initial conditions to form a self-consistent model de- 
scribed in section above. 

0.7 (1-96) Bugs fixed from previous version. 

0.9 (2-96) Completed version used for test case 3. Modified which data 
fields are dumped. Added max/min trace values. Added point trace 
values. 

0.9.1 (2-96) Restructured code. Put movie calls into routine do movie. 
Added restart capabilities, single SIZE parameter file. Added smart 
data directory. Released as solwnd_l . 0. 


13 User-friendly features 

solwnd has some user-friendly features making it a convenient code to work 
with. 


1. problem size file. The file SIZE_sw sets the problem size and only 
needs to be altered for different problem sizes. 

parameter (nx=132, ny=36, nz=36) 
c set nxmax to be the largest of (nx,ny,nz) 
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c 


set nymax to be the 2nd largest of (nx,ny,nz) 
parameter (nxmax=132, nymax=36) 


As indicated, nxmax and nymax need to be appropriately set. 
% input file. An input file is provided for the parameters. 


file solwnd_i .0. in 


/ disk6/ data/solwnd/ test/ 
sw_AAAA 
0 : if restart 

/disk6/data/ solwnd/test/sw_1005 . dat 


header (7 chars) 


1000 

200 

1 


nsteps : final iteration step 

idump : number of iterations between dump 

nstart : iteration starting counter 


0 . 25e+0 

crn 

: c our ant number 

1.00e-l 

bx_pert : delta 

B_0x, initial delta B in x direction 

0.10e+0 

eta__b 

: delta B rms/BO 

0 . 30e+0 

eta_perp: delta B_perp/B0 

0 . 30e+0 

eta_v 

: delta V_x rms 

l,80e+0 

fmach 

: Sonic Mach number 

1 . 30e+0 

fmacha 

: Alfven Mach number 

O.OOe+O 

xkx 

: B0 wavenumber in x direction 

1.00e+0 

xky 

: B0 wavenumber in y direction 

1.00e+0 

xkz 

: B0 wavenumber in z direction 

1.67e+0 

gamma 

: ratio of specific heats 

1.00e+0 

fk : shear layer thickness factor (1-2 is good) 

1.30e+l 

omega 

: frequency of incoming waves 

1.00e+l 

lambda 

: time constant for inflow ramp 

4.00e+0 

xlen 

: length of system in units of length 

1.00e+0 

ylen : 

length of system in units of length 

l.OOe+O 

zlen 


t 

lsynch 

: fct synchronization 

t 

lhydro 

: switch to turn on coupled hydro 

t 

lfield 

: switch to turn on magnetic field solver 


Line 1 is the directory in which all data, trace and movie files will be 
placed. The directory must be terminated with a 

Line 2 is the case name that will be prefixed to all the data, trace and 
movie files. A useful convention is sw_ tagging the AAAA which provide 
ample opportunity to vary parameters and successive runs can be la- 
beled AAAB, AAAC, AAAD etc. Data files (field files) are created with 
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the header and with 1000 (e.g. sw_AAAA1000 . dat , sw_AAAA1001 . dat , 
sw_AAAA1002.dat etc.) added to the dump number. Thus an Is or 
other unix sort lists them in order. 

Line 3 is a flag; if set to 1 indicates that the file listed in the next line 
must be read to initialize the data. The code begins the run with the 
step and time of the restart file. 

Line 4 is the restart file name. Note that it can be in a different 
directory from the present run. Point trace values are either created 
in a . trcl or . trc2 files or if they exist in the data directory they are 
appended to these files. 

The rest of the lines are either self explanatory or discussed in the text 
previously. 

3. data directory. As indicated above the data directory is written in 
with all data, trace and movie files. 

4. restart. As indicated above the code can restart from a previous run. 

5. movie files. The routine do .movie in the code creates subfields which 
are then passed to the routine dump .movie. New fields that are needed 
can be added to the ones already present or the present ones eliminated 
by editing do_movie. The program dump .movie . c is a code written by 
A. Malagoli (U. Chicago) for viewing the movie files. This program is 
invoked by, for instance, 

Movie -z 2 2 -p 5 movie_file_name 

The -p option indicates the color map. The -z option indicates a 
zoom of 2 in both directions. 

6. time step CPU counters . If -SGI is chosen in the compile line 
then a CPU time step counter is used. 

7. data, max/min and point trace values. Pull field data of density, 
pressure, 3 components of momentum, total energy and 3 components 
of magnetic fields are saved every ndump steps in .dat files. Every 
time step, time, max/min values of sonic and Alfvenic Mach number, 
density and pressure are saved in .trcl files. Point values of each of 
the primary variables of density, momenta, energy and magnetic fields 
are put in an unformatted file .trc2. The sample points are at the 
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stations, 


x — L x j 4? L x / 2) 3L x / 4, L x , 

located at y = L y / 2. For each of these stations 5 points in the vertical 
(z) direction at, 

z — L z j 4 2e, L z j 4 — e, L z /4^ / 4 4- e, L Z J 4 + 2e, 

are sampled, e is presently chosen as 4 grid points. 


14 Supplementary programs 

Several IDL (Interactive Data Language, Research Systems Inc.) routines 
have been written to analyze and display results from solwnd simulations. 

1. solwnd. pro reads the full field data. Depending on mflag it displays 
information about the max and min of sonic and Alfvenic Mach num- 
bers. It is invoked by, 

solwnd, if ile, mf lag, time , nx, ny ,nz,rho,p, vx,vy,vz,e,bx, by, bz,mach 

if ile is the file name (including path) of the solwnd field data files, 
e is the total energy, internal + kinetic + magnetic, mach is the three- 
dimensional sonic mach number field. Rest of the variables have their 
usual meaning. All the fields are of size (nx,ny,nz). 

2. solwnd_trc . pro reads the file case.trcl and displays the max/min of 
sonic and Alfvenic Mach numbers, density and pressure. It is invoked 
by, 

solwnd_trc,ifile,n,a 

if ile is the file name (including path) of the solwnd case.trcl file, 
n is number of time step data records (nsteps. a is of dimension 
(9,nsteps). 

a(0,*) is t 

a(l,*) is max M(x,y,z) 
a(2,*) is min M(x,y,z) 
a(3,*) is max MA(x,y,z) 
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a(4,*) is min MA{x,y,z) 
a(5,*) is max p(x,y,z) 
a(6 , *) is min p(x , y, z) 
a(7,*) is max p(x,y,z) 
a(8,*) is min p(x,y, z) 

Upon invocation, solwncLtrc . pro reads the trace file and then prompts 
the user for which variable is to be plotted in a composite max/min 
plot (where the plots for the Mach numbers also indicate the unity 
value) as follows, 

iplot => 0 skip plot 
iplot => 1 mach vs i 
iplot => 2 mach vs t 
iplot => 3 macha vs t 
iplot => 4 rho vs t 
iplot => 5 p vs t 


3. solwnd_trc2.pro reads the file case.trc2 (point data file) and plots 
the data vs. step or if “time” is available (from previous run of 
solwnd_trc . pro for instance) then plots the data vs. time. User 
is queried for what station (x location) and what height (z location) 
is plot of p,p, V x , V y , V z , B x , By, or B z vs. time is desired. The routine 
cycles through choices until an exit is indicated with a negative value. 

4. solwnd_arrow . pro reads the same data files as solwnd.pro but pro- 
duces various arrow plots. It is invoked as, 

solwnd_arrow , f ile , j cut , iplot , ires 

j cut is the y— plane in which the x—z view is desired, ires is 1 , 2, 4 
etc . indicating either full , half, or quarter resolution etc. The data is 
sub-sampled for drawing the arrow plots. For n x > 32 the arrow plots 
are too crowded to see them as arrows, although interesting effects can 
still be produced, iplot indicates one of the following, 

iplot => 1 (B x — Bq,B z ); (B x , 0) in x — z plane (2 plots/page) 
iplot => 2 (0, B z ); (0 ,B y ) in x — z plane (2 plots/page) 
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iplot => 3 (V x — Vo, V z ); (V x , 0) in x — z plane (2 plots/page) 
iplot => 4 (0, V z ); (0, V y ) in x — z plane (2 plots/page) 
iplot => 5 (Vy,V g ); (1^,0); (0,V Z ) in x - z plane (3 plots/page) 
iplot => 6 ( B y ,B z ); (B y , 0); (0 ,B Z ) in x — z plane (3 plots/page) 
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